#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <stdbool.h>

#include <inttypes.h>

uint32_t as_uint32(double f)
{
   union
   {
      double d;
      uint32_t i[2];
   } u = {f};
   return u.i[0];
}

int32_t as_int32(double f)
{
   union
   {
      double d;
      uint32_t i[2];
   } u = {f};
   return (int)(u.i[0]);
}

uint32_t convert_to_uint32(double x)
{
   double y = x + 0x1.8p52; // 0x1.8p52= 1.5 * 2^52
   return as_uint32(y);
}

int32_t convert_to_int32(double x)
{
   double y = x + 0x1.8p52; // 0x1.8p52= 1.5 * 2^52
   return as_int32(y);
}

// return x / 2^64, here x < 2^52
double div_2p64(uint64_t x)
{
   union
   {
      uint64_t i;
      double d;
   } u;
   u.i = ((1023lu - 12) << 52) | x;
   return (u.d - 0x1.0p-12);
}

void test_convert_to_uint32()
{
   double d;
   uint32_t n;
   printf("\nrun %s()\n", __func__);
   double ary[] =
       {
           1.0, 1.49, 1.5, 1.51, 2.0, 2.49, 2.5, 2.51,
           2147483648.0, 2147483648.49, 2147483648.50, 2147483648.51,
           4294967295.0, 4294967295.49, 4294967295.50, 4294967295.51};
   for (int i = 0; i < sizeof(ary) / sizeof(ary[0]); i++)
   {
      n = convert_to_uint32(ary[i]);
      printf("int(%f)=%u\n", ary[i], n);
   }
}

static double uint32_t_to_double(uint32_t n)
{
   union
   {
      uint32_t i[2];
      double d;
   } u;
   u.i[0] = n;
   u.i[1] = ((1023 + 52) << 20);
   return (u.d - 0x1.0p52);
}

double int64_to_double(int64_t i)
{
   return (double)(i);
}

/* -mavx
int64_to_double:
   vxorpd	xmm0, xmm0, xmm0
   vcvtsi2sdq	xmm0, xmm0, rdi
   ret
*/

double uint64_to_double(uint64_t i)
{
   return (double)(i);
}

/* __m128d _mm_cvtsi64_sd (__m128d a, __int64 b)
 Convert the signed 64-bit integer b to a double-precision (64-bit) floating-point element, //
 store the result in the lower element of dst, and copy the upper element from a to the
 upper element of dst.
*/

/* uint64_to_double 汇编代码

uint64_to_double:
   test	rdi, rdi
   js	.L14　　           ； 大于等于2^63转less0
   pxor	xmm0, xmm0　   ； xmm0=0
   cvtsi2sdq	xmm0, rdi
   ret
less0:
   mov	rax, rdi   ; rax=i
   pxor	xmm0, xmm0 ; xmm0=0
   shr	rax,1      ; rax= i 逻辑右移１位, rax= i/2;
   and	edi, 1　　　; 清edi 高63位
   or	rax, rdi      ; rax 最低位置0或者1;  rax= round(i/2)
   cvtsi2sdq	xmm0, rax ; xmm0= round(i/2)
   addsd	xmm0, xmm0      ; xmm0 = round(i/2)*2
   ret
*/

void test_uint64_to_double()
{
   printf("\nrun %s()\n", __func__);

   uint64_t m = ((1ull) << 53) - 1;
   double f = (double)m;
   printf("m=%" PRIu64 "=%lf\n", m, f);
   uint64_t i0 = (m << 2) + 0;
   uint64_t i1 = (m << 2) + 1;
   uint64_t i2 = (m << 2) + 2;
   uint64_t i3 = (m << 2) + 3;
   printf("m=%" PRIu64 "=%lf\n", i0, uint64_to_double(i0));
   printf("m=%" PRIu64 "=%lf\n", i1, uint64_to_double(i1));
   printf("m=%" PRIu64 "=%lf\n", i2, uint64_to_double(i2));
   printf("m=%" PRIu64 "=%lf\n", i3, uint64_to_double(i3));

   uint64_t i8 = (m << 3) + 4;
   printf("m=%" PRIu64 "=%lf\n", i8, uint64_to_double(i8));
}

void test_uint32_to_double()
{
   printf("\nrun %s()\n", __func__);

   uint32_t ns[] = {0, 1, 2, 4, 8, 16, 1024, 2048, 4095, 4096};

   for (size_t i = 0; i < sizeof(ns) / sizeof(ns[0]); i++)
   {
      uint32_t n = ns[i];
      double f=uint32_t_to_double(n);
      printf("%u=%lf\n", n, f);
   }
}

void test_convert_to_int32()
{
   printf("\nrun %s()\n", __func__);

   double d;
   int n;

   double ary[] = {
       -1.0, -1.49, -1.5, -1.51, -2.0, -2.49, -2.5, -2.51,
       1.0, 1.49, 1.5, 1.51, 2.0, 2.49, 2.5, 2.51};

   for (int i = 0; i < sizeof(ary) / sizeof(ary[0]); i++)
   {
      n = convert_to_uint32(ary[i]);
      printf("int(%f)=%d\n", ary[i], n);
   }
}

// double div_2p64(uint64_t x)
void test_div_2p64()
{
   printf("\nrun %s()\n", __func__);
   uint64_t ary[] = {0, 1, 2, 3, 4, 127, 128, 254, 255, 256, 4000, 4095};
   for (int i = 0; i < sizeof(ary) / sizeof(ary[0]); i++)
   {
      double r0 = (double)(ary[i]) * 0x1.0p-64;
      double r1 = div_2p64(ary[i]);
      printf("ref     =%" PRIu64 "/2^64=%.20lf=%la\n", ary[i], r0, r0);
      printf("div_2p64=%" PRIu64 "/2^64=%.20lf=%la\n\n", ary[i], r1, r1);
   }
}

void test_mod_128()
{
   printf("\nrun %s()\n", __func__);
   int min= -1022 * 128+1;
   int max=  1022 * 128-1;
   int mask=128-1;
   bool ok=true;
   int k0,k1,idx0,idx1;
   for (int i=min;i<=max;i++)
   {
      k0=(i>>7);
      idx0= i - k0*128;

      idx1= (i&mask);
      k1  = ((i-idx1)>>7 );
      if ( k0!=k1 || idx0!=idx1)
      {
         printf("i=%d, k0=%d,idx0=%d,k1=%d,idx1=%d\n",i,k0,idx0,k1,idx1);
         ok=false;
         break;
      }
   }
   if ( ok)
      printf("ok\n");
}


int main(int argc, char *argv[])
{
   test_mod_128();
   test_uint32_to_double();
   test_uint64_to_double();
   test_convert_to_uint32();
   test_convert_to_int32();
   test_div_2p64();
   return 0;
}
